/*                                                            02.Feb.2010  v2.7
  =============================================================================

                          U    U   GGG    SSSS  TTTTT
                          U    U  G       S       T
                          U    U  G  GG   SSSS    T
                          U    U  G   G       S   T
                           UUU     GG     SSS     T

                   ========================================
                    ITU-T - USER'S GROUP ON SOFTWARE TOOLS
                   ========================================


       =============================================================
       COPYRIGHT NOTE: This source code, and all of its derivations,
       is subject to the "ITU-T General Public License". Please have
       it  read  in    the  distribution  disk,   or  in  the  ITU-T
       Recommendation G.191 on "SOFTWARE TOOLS FOR SPEECH AND  AUDIO 
       CODING STANDARDS".
       =============================================================


MODULE:         EID.C, ERROR INSERTION DEVICE

ORIGINAL BY:    

                Rudolf Hofmann                                
                Advanced Development Digital Signal Processing
                PHILIPS KOMMUNIKATIONS INDUSTRIE AG           
                Kommunikationssysteme                         
                Thurn-und-Taxis-Strasse 14                    
                D-8500 Nuernberg 10 (Germany)                 
                                                              
                Phone : +49 911 526-2603                      
                FAX   : +49 911 526-3385                      
                EMail : hf@pkinbg.uucp

DESCRIPTION:      Error Insertion Device.
                  - Generates random or bursty error patterns
                    with bit error rates from 0.00% up to 50%
                  - disturbes data bits according these error patterns
                  - allows implementation of "frame erasure" with a
                    certain frame erasure rate (also random or bursty
                    nature selectable)

NOTES:         a) The format of I/O-data is explained in the header of
                  each function. There is also some information given 
                  concerning the algorithm of the function.

               b) Generation of bit error pattern and disturbance of the
                  individual bits have been splitted up into two functions.
                  This allows the load of error patterns (eg, from a file), 
                  which have been generated by other software
                  (i.e. with respect to certain channel models)

               c) The generation of the error patterns here is done with
                  the aid of a Gilber Elliot channel model, which allows the
                  production of random errors as well as bursty errors
                  by defining a simple 'correlation factor' (gamma).
                  For more details see description of function
                  "GEC_init" (see also separate ITU-T document).

FUNCTIONS       : - open_eid (double ber, double gamma)
                     Creates a struct containing arrays and variables
                     for the EID. Initializes the error pattern generator,
                     which produces error patterns defined by ber, gamma:
                     Input: ber   = bit error rate (0 =<  ber  < 0.5 )
                     Input: gamma = burst factor   (0 =< gamma < 0.99)
                                    the bursty nature of the error
                                    pattern increases with gamma
                                    (gamma = 0: equally distributed errors
                                     gamma = 1: totally bursty errors)

                  - BER_generator (SCD_EID *EID, long lseg, short *EPbuff)
                     Generates a bit error sequence with properties,
                     defined by "open_eid".
                     Input : EID    = pointer to EID-struct
                     Input : lseg   = number of error bits to be generated
                                      in one call.
                     Output: EPbuff = array, containing the error pattern

                  - BER_insertion (long lseg, short *ibuff,
                                              short *obuff, short *EPbuff)
                     Disturbes the input data bits according the error
                     pattern, stored in EPbuff[].
                     Input : lseg    = number of data bits in ibuff[]
                     Input : ibuff[] = array, containing "lseg" data bits
                                       (dimension: lseg+1;
                                        ibuff[0] = SYNC-word )
                     Input : EPbuff[]= array, containing the error pattern
                                       (dimension: lseg)
                     Output: obuff[] = array, containing the disturbed
                                       input data bits
                                       (dimension: lseg+1;
                                        obuff[0] = SYNC-word )
		  **	 
		  **  Include the following routine
		  **  G. Schroeder 04.08.1993
		  **	 
                  -  FER_generator_random(SCD_EID *EID)
                     Frame erasure module.
                     Input : EID     = pointer to EID-struct

                  - FER_module(SCD_EID *EID, long  lseg, short *ibuff,
                                                         short *obuff)
                     Frame erasure module.
                     Input : EID     = pointer to EID-struct
                     Input : lseg    = number of data bits in ibuff[]
                     Input : ibuff[] = array, containing "lseg" data bits
                                       (dimension: lseg+1;
                                        ibuff[0] = SYNC-word )
                     Output: obuff[] = array, containing the input data
                                       bits or all bits erased
                                       (dimension: lseg+1;
                                        obuff[0] = SYNC-word or
                                                   marked SYNC-word)
                  - close_eid (SCD_EID *EID)
                     Frees memmory allocated to the EID state variables' 
                     buffer. 

		   - open_burst_eid(long index);

		   - FER_generator_burst(BURST_EID *state);
                   - reset_burst_eid(BURST_EID *burst_eid);

 HISTORY:
  28.Feb.92 v1.0 1st UGST version
  20.Apr.92 v2.0 Modifications on the RNG
  24.Apr.92 v2.1 Modifications for portability issues and inclusion of
                 the PORT_TEST compilation clausule.
  24.Apr.92 v2.2 Inclusion of close_eid.
  04.Aug.93 v2.3 Inclusion of burst frame erasure routines by 
                 <gerhard.schroeder@ties.itu.ch>
  06.Jun.95 v2.4 Changed self-documentation of module BER_insertion to 
                 align the softbit definitions of the EID module 
                 with the EID hardware and Rec.G.192. Code not
                 changed <simao@ctd.comsat.com>
  15.Aug.97 v2.5 Added a function to reset the Belcore mode state
                 variable, reset_burst_eid() <simao>
  12.Jan.07 v2.6 open_burst_eid modified and function bfer_comp added 
                 to extend Bellcore burst model resolution and operating
                 range to [0.5-30%]. <J.Sv. Ericsson> 
  02.Feb.10 v2.7 Modified maximum string lenght for filenames (y.hiwasaki)
  =============================================================================
*/

 
 
/* ********************* BEGIN OF EXAMPLES ON USAGE *********************** */
 
/* ==========================================================================
   Below two example programs are presented, how to use the EID-functions,
   one for insertion of bit errors and the other one for frame erasure.

   VERY IMPORTANT NOTE: If the internal EID-states are not saved to file at 
      the end of the program, identical error patterns will be produced when 
      the program is re-started.
      Since file-I/O must be implemented by the user, examples for this
      are given in the EIDDEMO.C demo program.

 ============================================================================
*/

/*============================================================================*/
/*                   Example Program for BIT ERROR INSERTION                  */
/*============================================================================*/
#ifdef A_Z_R_J_B_L_X                    /* to avoid compilation of examples */

main(int argc, char** argv)
{
   SCD_EID *BEReid;                     /* pointer to EID-structure */
   char ifile[MAX_STRLEN],ofile[MAX_STRLEN]; /* input/output file names */
   FILE *ifilptr,*ofilptr;              /* input/output file pointer */
   int    EOF_detected;
 
   double ber,gamma;                    /* bit error rate, burst factor */
   long   lseg,lread;                   /* length of one transmitted frame */
   double ber1;                         /* return values from BER_generator */
   double dstbits;                      /* count total distorted bits */
   double prcbits;                      /* count number of processed bits */

   short *xbuff,*ybuff,shvar;           /* pointer to bit-buffer */
   short *EPbuff;                       /* pointer to error pattern-buffer */
   short SYNCword,i;
 
    clock_t t1,t2;
    double t;
 
/* ------------------------------- Statements ------------------------------- */
 
   .... Put your statements here ....
 
  GET_PAR_S(1,"_File with input bit stream: ............... ",ifile);
   ifilptr = fopen(ifile,RB);           /* open input file */
   if (ifilptr==NULL)
      HARAKIRI("    Could not open input file",1);
 
        /* ------------------------------------------------------------------ */
        /* Find SYNC-word and frame length on file */
        /* ------------------------------------------------------------------ */
        lread = fread(&SYNCword,2,1,ifilptr);
        if (SYNCword!=(short)0x6B21)
          HARAKIRI("    First word on input file not the SYNC-word (0x6B21)",1);
        lseg = 0;
        i    = 0;
        while(lseg==0)
        {
            lread = fread(&shvar,2,1,ifilptr);
            if (lread==0)
               HARAKIRI("    No next SYNC-word found on input file",1);
            if (shvar==SYNCword)
                lseg=i;
            i++;
        }
        fseek(ifilptr,0L,0);            /* filepointer back */
 
 
  GET_PAR_S(2,"_File for disturbed bitstream: ............. ",ofile);
   ofilptr = fopen(ofile,WB);
   if (ofilptr==NULL)
      HARAKIRI("    Could not create output file",1);
 
 
/* -------------------------------------------------------------------------- */
/*                       Initialize EID (bit error rate) */
/* -------------------------------------------------------------------------- */
  GET_PAR_D(3,"_bit error rate     (0.0 ... 0.50): ........ ",ber);
  GET_PAR_D(4,"_burst factor       (0.0 ... 0.99): ........ ",gamma);
 
     BEReid = open_eid(ber,gamma);      /* Open EID */
     if (BEReid==(SCD_EID *)0)
        HARAKIRI("    Could not create EID for bit errors!?",1);
 
 
/* -------------------------------------------------------------------------- */
/*                       Allocate memory for I/O-buffer */
/* -------------------------------------------------------------------------- */
  printf("_found frame length on input file: ......... %d\n",lseg);
 
/* -------------------------------------------------------------------------- */
/* a) buffer for input bit stream: */
    xbuff = (short *)malloc((1+lseg)*sizeof(short));
    if (xbuff==(short *)0)
       HARAKIRI("    Could not allocate memory for input bit stream buffer",1);
 
/* -------------------------------------------------------------------------- */
/* b) buffer for output bit stream: */
    ybuff = (short *)malloc((1+lseg)*sizeof(short));
    if (ybuff==(short *)0)
       HARAKIRI("    Could not allocate memory for output bit stream buffer",1);
 
/* -------------------------------------------------------------------------- */
/* c) buffer for storing the error pattern: */
    EPbuff = (short *)malloc((lseg)*sizeof(short));
    if (EPbuff==(short *)0)
       HARAKIRI("    Could not allocate memory for error pattern buffer",1);
 
/* -------------------------------------------------------------------------- */
/*                           Now process input file */
/* -------------------------------------------------------------------------- */
/* initialize counters fro computing bit error rate/frame erasure rate */
    dstbits = 0.0;                      /* number of distorted bits */
    prcbits = 0.0;                      /* number of processed bits */

    /* ---------------------------------------------------------------------- */
    /*      Read file with input bit stream (segments of lenght "lseg+1") */
    /* ---------------------------------------------------------------------- */
    EOF_detected=0;                     /* Flag, to mark END OF FILE */
    t1 = clock();                       /* measure CPU-time */
 
    while( (lread=fread(xbuff,2,lseg+1,ifilptr)) == lseg+1)
    {
      if (xbuff[0]==SYNCword && EOF_detected==0)
      {
        /* ------------------------------------------------------------------ */
        /*                       Generate Error Pattern */
        /* ------------------------------------------------------------------ */
        ber1 = BER_generator(BEReid, lseg,EPbuff);
          dstbits += ber1;              /* count number of disturbed bits */
          prcbits += (double)lseg;      /* count number of processed bits */

        /* ------------------------------------------------------------------ */
        /*     Modify input bit stream according the stored error pattern */
        /* ------------------------------------------------------------------ */
        BER_insertion(lseg+1,xbuff,ybuff,EPbuff);
 
 
        /* ------------------------------------------------------------------ */
        /*               Write disturbed bits to output file */
        /* ------------------------------------------------------------------ */
        lread=fwrite(ybuff,2,lseg+1,ofilptr);
      }
 
      else                              /* if the next SYNC-word is missed */
      {
         EOF_detected=1;
      }
    }
 
    if (EOF_detected==1)
        printf("   --- end of file detected (no SYNCword match) ---\n");
    printf("\n");
 
    t2 = clock();
    t  = (t2 - t1) / (double) CLOCKS_PER_SEC;
    printf("CPU-time %f sec\r\n\n",t);
 
/* -------------------------------------------------------------------------- */
/*                 Print message with measured bit error rate */
/* -------------------------------------------------------------------------- */
    if (prcbits>0)
    {
        printf("measured bit error rate    : %f  (%ld of %ld bits distorted)\n",
                dstbits/prcbits,(long)dstbits,(long)prcbits);
    }
 
    exit(0);
}
 
 
 
 
 
 
 
 
 
/* ========================================================================== */
/*                      Example program for FRAME ERASURE                     */
/* ========================================================================== */
/* -------------------------------------------------------------------------- */
/*                                Main Program                                */
/* -------------------------------------------------------------------------- */

main(int argc, char** argv)
{
   SCD_EID *FEReid;                     /* pointer to EID-structure */
   char ifile[MAX_STRLEN],ofile[MAX_STRLEN]; /* input/output file names */
   FILE *ifilptr,*ofilptr;              /* input/output file pointer */
   int    EOF_detected;
 
   double fer,gamma;                    /* frame erasure rate; burst factor */
 
   long   lseg,lread;                   /* length of one transmitted frame */
   double fer1;                         /* return value from FER_insertion */
   double ersfrms;                      /* total distorted frames */
   double prcfrms;                      /* number of processed frames */
   short *xbuff,*ybuff,shvar;           /* pointer to bit-buffer */
   short SYNCword,i;
 
   clock_t t1,t2;
   double t;
 
/* ------------------------------- Statements ------------------------------- */
 
  .... Put your statements here ....
 
  GET_PAR_S(1,"_File with input bit stream: ............... ",ifile);
   ifilptr = fopen(ifile,RB);
   if (ifilptr==NULL)
      HARAKIRI("    Could not open input file",1);
        /* ------------------------------------------------------------------ */
        /* Find SYNC-word and frame length on file */
        /* ------------------------------------------------------------------ */
        lread = fread(&SYNCword,2,1,ifilptr);
        if (SYNCword!=(short)0x6B21)
          HARAKIRI("    First word on input file not the SYNC-word (0x6B21)",1);
        lseg = 0;
        i    = 0;
        while(lseg==0)
        {
            lread = fread(&shvar,2,1,ifilptr);
            if (lread==0)
               HARAKIRI("    No next SYNC-word found on input file",1);
            if (shvar==SYNCword)
                lseg=i;
            i++;
        }
        fseek(ifilptr,0L,0);            /* filepointer back */
 
 
  GET_PAR_S(2,"_File for disturbed bitstream: ............. ",ofile);
   ofilptr = fopen(ofile,WB);
   if (ofilptr==NULL)
      HARAKIRI("    Could not create output file",1);
 
  GET_PAR_D(3,"_frame erasure rate (0.0 ... 0.50): ........ ",fer);
  GET_PAR_D(4,"_burst factor       (0.0 ... 0.99): ........ ",gamma);
 
/* -------------------------------------------------------------------------- */
/*                               Initialize EID */
/* -------------------------------------------------------------------------- */
     FEReid = open_eid(fer,gamma);      /* Open EID for frame erasure */
     if (FEReid==(SCD_EID *)0)
        HARAKIRI("    Could not create EID for frame erasure module!?",1);
 
 
/* -------------------------------------------------------------------------- */
/*                       Allocate memory for I/O-buffer */
/* -------------------------------------------------------------------------- */
  printf("_found frame length on input file: ......... %d\n",lseg);
 
/* -------------------------------------------------------------------------- */
/* a) buffer for data from input file: */
    xbuff = (short *)malloc((1+lseg)*sizeof(short));
    if (xbuff==(short *)0)
       HARAKIRI("    Could not allocate memory for input bit stream buffer",1);
 
/* -------------------------------------------------------------------------- */
/* b) buffer for output bit stream: */
    ybuff = (short *)malloc((1+lseg)*sizeof(short));
    if (ybuff==(short *)0)
       HARAKIRI("    Could not allocate memory for output bit stream buffer",1);
 
/* -------------------------------------------------------------------------- */
/*                           Now process input file */
/* -------------------------------------------------------------------------- */
/* initialize counters fro computing bit error rate/frame erasure rate */
    ersfrms = 0.0;                      /* number of erased frames */
    prcfrms = 0.0;                      /* number of processed frames */

    /* ---------------------------------------------------------------------- */
    /*      Read file with input bit stream (segments of lenght "lseg+1") */
    /* ---------------------------------------------------------------------- */
    EOF_detected=0;                     /* Flag, to mark END OF FILE */
    t1 = clock();                       /* measure CPU-time */
    while( (lread=fread(xbuff,2,lseg+1,ifilptr)) == lseg+1)
    {
      if (xbuff[0]==SYNCword && EOF_detected==0)
      {
        /* ------------------------------------------------------------------ */
        /*                            Frame erasure */
        /* ------------------------------------------------------------------ */
        fer1 = FER_module(FEReid,lseg+1,xbuff,ybuff);
          ersfrms += fer1;              /* count number of erased frames */
          prcfrms += (double)1;         /* count number of processed frames */
 
        /* ------------------------------------------------------------------ */
        /*                    Write (erased) frames to file */
        /* ------------------------------------------------------------------ */
        lread=fwrite(ybuff[0],2,lseg+1,ofilptr);
      }
 
      else                              /* if the next SYNC-word is missed */
      {
         EOF_detected=1;
      }
    }
 
    if (EOF_detected==1)
        printf("   --- end of file detected (no SYNCword match) ---\n");
    printf("\n");
 
    t2 = clock();
    t  = (t2 - t1) / (double) CLOCKS_PER_SEC;
    printf("CPU-time %f sec\r\n\n",t);
 
/* -------------------------------------------------------------------------- */
/*                 Print message with measured bit error rate */
/* -------------------------------------------------------------------------- */
    if (prcfrms>0)
    {
        printf("measured frame erasure rate: %f  (%ld of %ld frames erased)\n",
                ersfrms/prcfrms,(long)ersfrms,(long)prcfrms);
    }
 
    exit(0);
}
#endif
/* *********************** END OF EXAMPLES ON USAGE *********************** */
 

#ifdef PORT_TEST
    int PORTABILITY_TEST_OPERATION = 1;
#endif

 
 
/* ......... Include of general definitions .........*/ 
#include <time.h>
#include <math.h>
#include <stdlib.h>
#include <ctype.h>
#include <stdio.h>

/* ......... Include of EID prototypes and definitions .........*/ 
#include "eid.h"

/* Local function prototypes and definitions .........*/ 
double EID_random ARGS((unsigned long *seed));
void update_EID_random ARGS((long len_register, long *shift_register));
long GEC_init ARGS((SCD_EID *EID, double  ber, double  gamma));
double bfer_comp(long index);
/* ......... Global variable: Bellcore Model Transition probability vector ......... */
/* ......... prob[0]=P0 updated from command line, using function bfer_comp......... */
double prob[MODEL_SIZE] = { 0.0023, 0.85, 0.825, 0.8, 0.775, 0.75, 0.725, 0.7, 0.6, 0.45, 0.0}; 

/*
 * ...................... BEGIN OF FUNCTIONS .........................
 */
 

/*
  ============================================================================ 

        SCD_EID *open_eid (double ber, double gamma);
        ~~~~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Allocate memory for EID struct and setup the transission matrix
        according the select bit error rate.

        Parameters:  
        ~~~~~~~~~~~
        ber: ..... (In) bit error rate 
        gamma: ... (In) burst factor 

        Return value:
        ~~~~~~~~~~~~~
        Returns a pointer to struct SCD_EID;

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        28.Feb.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================
*/
SCD_EID *open_eid (ber, gamma)
double ber,gamma;
{
    SCD_EID *EID;
    time_t  t1;
 

  /* Allocate EID structure */
   if ((EID=(SCD_EID *)malloc(sizeof(SCD_EID)))==0L) return((SCD_EID *)0);

  /* Preset of the random generator seed with current system time */
#ifdef PORT_TEST
    t1 = 314159265;
#else
    time(&t1);
#endif
    EID->seed = (unsigned long) t1;
 

  /* Initialize Gilbert-Elliot Channel model */
    if (GEC_init(EID,ber,gamma)==0L) return((SCD_EID *)0);
 
  /* Store ber/gamma in struct */
        EID->usrber   = ber;            /* user defined bit error rate */
        EID->usrgamma = gamma;          /* user defined correlation factor */

  /* Return initialized structure */
        return(EID);
}  
/* ....................... End of open_eid() ....................... */
 
 

/*
  ============================================================================ 

        void close_eid (SCD_EID *EID);
        ~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Release the memory previously allocated by open_eid() for 
        an EID struct.

        Parameters:  
        ~~~~~~~~~~~
        eid: ..... (InOut) pointer to struct SCD_EID 

        Return value:
        ~~~~~~~~~~~~~
        None.

        Author: <tdsimao@venus.cpqd.ansp.br>
        ~~~~~~~

        History:
        ~~~~~~~~
        23.Apr.92 v1.0 Release of 1st version <tdsimao@venus.cpqd.ansp.br>

 ============================================================================
*/
void close_eid (EID)
    SCD_EID *EID;
{
   long i;


   /* Free state transition matrix columns */
    for (i=0; i<EID->nstates; i++)
       free((char *)EID->matrix[i]);

   /* Free state transition matrix rows */
    free((char *)EID->matrix);

  /* Free memory of bit error array */
    free((char *)EID->ber);

   /* Free EID structure */
   free((char *)EID);
}  
/* ....................... End of close_eid() ....................... */
 
 

/*                   
  ============================================================================ 

        double BER_generator (SCD_EID *EID, long lseg, short *EPbuff);
        ~~~~~~~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Generates bit error pattern according to the selected channel
        model.

        Parameters:  
        ~~~~~~~~~~~
        EID: ...... (In/Out) struct with channel model 
        lseg: ..... (In)     length of current frame 
        EPbuff: ... (Out)    bit error pattern (softbits) 

        Return value:
        ~~~~~~~~~~~~~
        Returns the bit error rate in the current frame as a double.

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        28.Feb.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================
*/
double BER_generator (EID, lseg,EPbuff)
SCD_EID *EID;
long    lseg;
short   *EPbuff;
{
   long i,n;                        /* value of random generator */
   double RAN, ber;                 /* aux. for computing bit error rate */
 

   /* Return if no samples are to be processed */
    if (lseg==(long)0) return(0.0);
 
   /* Generate random bits */
    ber = 0.0;
    for (i=0; i<lseg; i++)
    {
        /* Get next random number */
        RAN = EID_random(&(EID->seed));
 
        /* ... check the new channel state */
        for (n=0; n<EID->nstates; n++)
        {
            if (RAN < EID->matrix[EID->current_state][n])
            {
                EID->current_state = n; /* goto to the selected state */
                n = EID->nstates;       /* -> aborts loop */
            }
        }
 

  /* 
   * ......... COMPUTE BIT ERROR IN CURRENT STATE .........
   */

        /* Get next random number */
        RAN = EID_random(&(EID->seed));
 
        /* If random number is below the threshold (bit error rate in 
           current state), insert soft decision information for 'error' */
        if (RAN < EID->ber[EID->current_state])
        {   EPbuff[i] = (short)0x0081;
            ber += 1.0;                 /* increment number of errors */
        }
        /* otherwise insert soft decision information for 'no error' */
        else
        {   EPbuff[i] = (short)0x007F;
        }
    }
    return(ber);                        /* return number of error bits */
} 
/* ....................... End of BER_generator() ....................... */

 
/*                                                                            
  ============================================================================ 

        void BER_insertion_stl92 (long lseg, short *xbuff, short *ybuff, 
        ~~~~~~~~~~~~~~~~~~~~~~~~  short *error_pattern);

        Description:
        ~~~~~~~~~~~~

        Disturbing an input bitstream according stored error patterns.

        ---------------------------------------------------------------

                         NOTES on Data Representation:
                         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        Each individual bit is stored in one location of a  short-array,
        where in location 0 the SYNC-word is  stored (SYNC-word =
        0x6B21, 0x6B22, ... 0x6B2F).  To allow the processing of
        'soft-bits' (see separate  ITU-T-paper) in future applications,
        the following  definitions have been made: 

        a) input signal: 'hard' bit '0' must be represented as 0x007F 
                         'hard' bit '1' must be represented as 0x0081 

        b) error pattern: 
        * probability for undisturbed transmission: 
            0x0001 ... 0x007F       (0x0001: lowest probability that the 
            or   1 ... 127 dec.      transmission has been indeed undisturbed)

        * probability for disturbed transmission: 
            0x00FF ... 0x0081       (0x00FF: lowest probability that the
            or 255 ... 129 dec.      transmission has been really disturbed)


        c) output signal computation: 

        For input '1' (0x0081): 
          * if error pattern 0x00FF ... 0x0081 (255 ... 129), then
                   output = 0x0001 ... 0x007F  (  1 ... 127);
          * if error pattern 0x0001 ... 0x007F (  1 ... 127), then 
                   output = 0x00FF ... 0x0081  (255 ... 129).

        For input '0' (0x007F): 
          * if error pattern 0x00FF ... 0x0081 (255 ... 129), then 
                   output = 0x00FF ... 0x0081  (255 ... 129);
          * if error pattern 0x0001 ... 0x007F (  1 ... 127), then
                   output = 0x0001 ... 0x007F  (  1 ... 127).
        ---------------------------------------------------------------


        Parameters:  
        ~~~~~~~~~~~

        lseg: ..... length of current frame (with SYNC-word)

        xbuff: .... buffer with input bitstream (length = "lseg")

        ybuff: .... buffer with output bitstream (length = "lseg")

        EPbuff: ... buffer with error pattern (without SYNC-word)
                    (length = "lseg-1") 

        Return value:
        ~~~~~~~~~~~~~
        None.

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        28.Feb.92 v1.0 Release of 1st version <hf@pkinbg.uucp>
        06.Jun.95 v1.1 Changed self-documentation to align the EID module 
                       with the EID hardware and Rec.G.192. Code not
                       changed <simao@ctd.comsat.com>
 ============================================================================
*/
void BER_insertion_stl92(lseg, xbuff,ybuff, EPbuff)
    long lseg;
    short *xbuff; 
    short *ybuff; 
    short *EPbuff;
{
   long i;
 
  /* FIRST STEP: Copy SYNC-word -> output buffer */
    ybuff[0] = xbuff[0];

  /* THEN: Process data bits */
    /*    for (i=1; i<lseg+1; i++) This was in the STL92, but this is wrong!
     */
    for (i=1; i<lseg; i++)
    {   
      if (xbuff[i]==(short)0x007F)
        {
          /* input bit = '0' ('hard bit' = 7F) */
          ybuff[i] = (short)(0x00FF) &   EPbuff[i-1];
        }
      else
        {  
          /* input bit = '1' ('hard bit' = 81) */
          ybuff[i] = (short)(0x00FF) & (~EPbuff[i-1] + 1);
        }
    }
}
/* ..................... End of BER_insertion_stl92() ..................... */

 
 
/*                                                                            
  ============================================================================ 

        void BER_insertion_stl96 (long lseg, short *xbuff, short *ybuff, 
        ~~~~~~~~~~~~~~~~~~~~~~~~  short *error_pattern);

        Description:
        ~~~~~~~~~~~~

        Disturbing an input bitstream according stored error
        patterns. The input (undisturbed) and output (disturbed)
        buffers have samples conforming to the bitstream
        representation description in Annex B of G.192. This bitstream
        has a synchronism (sync) header at frame boundaries, which
        consists of a synchronism (SYNC) word followed by a 16-bit,
        2-complement word with the number of softbits in the frame
        (i.e., the frame length). The error-pattern does not have the
        sync header.

        This function is aliased to BER_insertion by default in
        eid.h. If the pre-STL96 version is desired, the module needs
        to be compiled with the symbol STL92 defined.

        ---------------------------------------------------------------

                         NOTES on Data Representation:
                         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        Each individual bit is stored in one location of a  short-array,
        where in location 0 the SYNC-word is  stored (SYNC-word =
        0x6B21, 0x6B22, ... 0x6B2F).  To allow the processing of
        'soft-bits' (see separate  ITU-T-paper) in future applications,
        the following  definitions have been made: 

        a) input signal: 'hard' bit '0' must be represented as 0x007F 
                         'hard' bit '1' must be represented as 0x0081 

        b) error pattern: 
        * probability for undisturbed transmission: 
            0x0001 ... 0x007F       (0x0001: lowest probability that the 
            or   1 ... 127 dec.      transmission has been indeed undisturbed)

        * probability for disturbed transmission: 
            0x00FF ... 0x0081       (0x00FF: lowest probability that the
            or 255 ... 129 dec.      transmission has been really disturbed)


        c) output signal computation: 

        For input '1' (0x0081): 
          * if error pattern 0x00FF ... 0x0081 (255 ... 129), then
                   output = 0x0001 ... 0x007F  (  1 ... 127);
          * if error pattern 0x0001 ... 0x007F (  1 ... 127), then 
                   output = 0x00FF ... 0x0081  (255 ... 129).

        For input '0' (0x007F): 
          * if error pattern 0x00FF ... 0x0081 (255 ... 129), then 
                   output = 0x00FF ... 0x0081  (255 ... 129);
          * if error pattern 0x0001 ... 0x007F (  1 ... 127), then
                   output = 0x0001 ... 0x007F  (  1 ... 127).
        ---------------------------------------------------------------


        Parameters:  
        ~~~~~~~~~~~

        lseg: ..... length of current frame (WITH sync header). This 
                    has been MODIFIED from the definition in the
                    pre-STL96 version of this function. Therefore lseg
                    is the number of softbits in the input bitstream
                    plus two (to account for the sync word and for the
                    length word).

        xbuff: .... buffer with input G.192 bitstream (length = "lseg")

        ybuff: .... buffer with output G.192 bitstream (length = "lseg")

        EPbuff: ... buffer with error pattern (without sync header)
                    (length = "lseg-2") 

        Return value:
        ~~~~~~~~~~~~~
        None.

        Author:
        ~~~~~~~
        Simao Ferraz de Campos Neto
        Comsat Laboratories                  Tel:    +1-301-428-4516
        22300 Comsat Drive                   Fax:    +1-301-428-9287
        Clarksburg MD 20871 - USA            E-mail: simao@ctd.comsat.com

        History:
        ~~~~~~~~
        07.Mar.96 v1.0 Created based on the STL92 version of this function, 
                       which did not comply with the bitstream data
                       representation in Annex B of G.192.
                       <simao@ctd.comsat.com>
 ============================================================================
*/
void BER_insertion_stl96(lseg, xbuff,ybuff, EPbuff)
    long lseg;
    short *xbuff; 
    short *ybuff; 
    short *EPbuff;
{
  register long i;
 
  /* FIRST STEP: Copy sync header to the output buffer */
  ybuff[0] = xbuff[0];
  ybuff[1] = xbuff[1];

  /* THEN: Process data bits */
  for (i=2; i<lseg; i++)
  {   
    if (xbuff[i]==(short)0x007F)
    {
      /* input bit = '0' ('hard bit' = 7F) */
      ybuff[i] = (short)(0x00FF) &   EPbuff[i-2];
    }
    else
    {  
      /* input bit = '1' ('hard bit' = 81) */
      ybuff[i] = (short)(0x00FF) & (~EPbuff[i-2] + 1);
    }
  }
}
/* ..................... End of BER_insertion_stl96() ..................... */

 
 
/*
  ============================================================================ 

        double FER_generator_random (SCD_EID *EID);

        Description:
        ~~~~~~~~~~~~

        Implementation of the `Frame Erasure' function:
        - computes the 'frame erasure pattern' 


        Parameters:   
        ~~~~~~~~~~~ 

        EID: ..... pointer to EID struct.

        Return value:
        ~~~~~~~~~~~~~
        Returns (double)1 if the current frame has been erased, and
        (double)0 otherwise. 

        Author: <hf@pkinbg.uucp>
        ~~~~~~~
	Modified by G. Schroeder 

        History:
        ~~~~~~~~
        28.Feb.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================ 
*/
double FER_generator_random(EID)
SCD_EID *EID; 
{
   long n;                        /* value of random generator */
   double RAN,fer;
 

   /* Get next random number */
   RAN = EID_random(&(EID->seed));
 
   /* See if another channel state has to be entered */
   for (n=0; n<EID->nstates; n++)
   {
     if (RAN < EID->matrix[EID->current_state][n])
     {
       EID->current_state = n; /* go to the selected state */
       n = EID->nstates;       /* -> aborts loop */
     }
   }
 
/*
   * ......... Compute bit error in current state .........
   */

   /* Get next random number */
   RAN = EID_random(&(EID->seed));
 
   /* if random number below the threshold (bit error rate in current state), 
      erase current frame totally */
   if (RAN < EID->ber[EID->current_state])
   {
     fer = 1.0;
   }
   else
   {
     fer = 0.0;
   }

  /* Return the frame erasure flag */
    return(fer);
}
/* ....................... End of FER_generator_random() ................. */


/*
  ============================================================================ 

        double FER_module_stl92 (SCD_EID *EID, long lseg, short *xbuff, 
        ~~~~~~~~~~~~~~~~~~~~~~~  short *ybuff);

        Description:
        ~~~~~~~~~~~~

        Implementation of the `Frame Erasure' function:
        - computes the 'frame erasure pattern' 
        - erases all bits in one frame according the current 
          state of the pattern generator. 

        The input buffer contains 'lseg' words, containing the data bits,
        preceeded by a SYNC-word.  If the frame should be erased (depending 
        on the frame  erasure pattern), all bits are set to 0x0000.
        Since a '0' is represented by 0x7F and a '1' is  represented 
        by 0x81, this corresponds to a total uncertainty about the true bit
        values. 

        In addition the lower 4 bits of the SYNC-word in the sync
        header are set to '0'. This makes it easier for the
        succeeding software to detect an erased frame.


        Parameters:   
        ~~~~~~~~~~~ 

        EID: ..... pointer to EID struct.
        lseg: .... length of current frame (including SYNC-word).
        xbuff: ... buffer with input bit stream xbuff[0](= SYNC-word) is
                   not processed.
        ybuff: ... buffer with output bit stream 

        Return value:
        ~~~~~~~~~~~~~
        Returns (double)1 if the current frame has been erased, and
        (double)0 otherwise. 

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        28.Feb.92 v1.0 Release of 1st version <hf@pkinbg.uucp>
        06.Jun.95 v1.1 Changed self-documentation to align the EID module 
                       with the EID hardware and Rec.G.192. Code not
                       changed <simao@ctd.comsat.com>

 ============================================================================ 
*/
double FER_module_stl92(EID, lseg,xbuff,ybuff)
SCD_EID *EID; 
long  lseg; 
short *xbuff; 
short *ybuff; 
{
   long i,n;                        /* value of random generator */
   double RAN,fer;
 

        /* Get next random number */
        RAN = EID_random(&(EID->seed));
 
        /* See if another channel state has to be entered */
        for (n=0; n<EID->nstates; n++)
        {
            if (RAN < EID->matrix[EID->current_state][n])
            {
                EID->current_state = n; /* go to the selected state */
                n = EID->nstates;       /* -> aborts loop */
            }
        }
 
/*
   * ......... Compute bit error in current state .........
   */

        /* Get next random number */
        RAN = EID_random(&(EID->seed));
 
        /* if random number below the threshold (bit error rate in current state), 
           erase current frame totally */
        if (RAN < EID->ber[EID->current_state])
        {
            ybuff[0] = xbuff[0] & 0x0000FFF0; /* modify SYNC-word */
            /* set all bits in current frame to 0 (which means 'not good - 
               not bad') */
            for (i=1; i<lseg; i++) ybuff[i] = (short)0;    
            fer = 1.0;
        }
        else
        {
        /* otherwise copy current frame into output buffer */
            ybuff[0] = xbuff[0];        /* store SYNC-word */
            for (i=1; i<lseg; i++)      /* store data */
                ybuff[i] = (short)(0x00FF) & (xbuff[i]);
            fer = 0.0;
        }

  /* Return the frame erasure flag */
    return(fer);
}
/* ....................... End of FER_module_stl92() ....................... */
 

/*
  ============================================================================ 

        double FER_module_stl96 (SCD_EID *EID, long lseg, short *xbuff, 
        ~~~~~~~~~~~~~~~~~~~~~~~  short *ybuff);

        Description:
        ~~~~~~~~~~~~

        Implementation of the `Frame Erasure' function:
        - computes the 'frame erasure pattern' 
        - erases all bits in one frame according the current 
          state of the pattern generator. 

        The input (undisturbed) and output (disturbed)
        buffers have samples conforming to the bitstream
        representation description in Annex B of G.192. This bitstream
        has a synchronism (sync) header at frame boundaries, which
        consists of a synchronism (SYNC) word followed by a 16-bit,
        2-complement word with the number of softbits in the frame
        (i.e., the frame length). Should the frame be erased (depending 
        on the frame  erasure pattern), all softbits are set to 0x0000.
        Since a '0' is represented by 0x7F and a '1' is  represented 
        by 0x81, this corresponds to a total uncertainty about the true bit
        values. 

        In addition the lower 4 bits of the SYNC-word in the sync
        header are set to '0'. This makes it easier for the succeeding
        software to detect an erased frame. The frame length word is
        copied unmodified to the output buffer.

        This function is aliased to FER_module by default in
        eid.h. If the pre-STL96 version is desired, the module needs
        to be compiled with the symbol STL92 defined.

        Parameters:   
        ~~~~~~~~~~~ 

        EID: ..... pointer to EID struct.
        lseg: .... length of current frame (including SYNC header).
        xbuff: ... buffer with undisturbed input G.192-compliant 
                   soft bitstream (length= lseg)
        ybuff: ... buffer with disturbed output G.192-compliant 
                   soft bitstream (length= lseg)

        Return value:
        ~~~~~~~~~~~~~
        Returns (double)1 if the current frame has been erased, and
        (double)0 otherwise. 

        Author:
        ~~~~~~~
        Simao Ferraz de Campos Neto
        Comsat Laboratories                  Tel:    +1-301-428-4516
        22300 Comsat Drive                   Fax:    +1-301-428-9287
        Clarksburg MD 20871 - USA            E-mail: simao@ctd.comsat.com

        History:
        ~~~~~~~~
        07.Mar.96 v1.0 Created based on the STL92 version of this function
                       by <hf@pkinbg.uucp>, which did not comply with
                       the bitstream data representation in Annex B of
                       G.192. <simao@ctd.comsat.com>

 ============================================================================ 
*/
double FER_module_stl96(EID, lseg,xbuff,ybuff)
SCD_EID *EID; 
long  lseg; 
short *xbuff; 
short *ybuff; 
{
  long i,n;                        /* value of random generator */
  double RAN,fer;
 

  /* Get next random number */
  RAN = EID_random(&(EID->seed));
 
  /* See if another channel state has to be entered */
  for (n=0; n<EID->nstates; n++)
  {
    if (RAN < EID->matrix[EID->current_state][n])
    {
      EID->current_state = n; /* go to the selected state */
      n = EID->nstates;       /* -> aborts loop */
    }
  }
 
  /*
   * ......... Compute bit error in current state .........
   */

  /* Get next random number */
  RAN = EID_random(&(EID->seed));
 
  /* if random number below the threshold (bit error rate in current state), 
     erase current frame totally */
  if (RAN < EID->ber[EID->current_state])
  {
    /* Modify SYNC-word and copy frame length word */
    ybuff[0] = xbuff[0] & 0x0000FFF0; 
    ybuff[1] = xbuff[1];

    /* set all bits in current frame to 0x0000 (which means 'not good - 
       not bad') */
    for (i=2; i<lseg; i++)
      ybuff[i] = (short)0x0000;    
    fer = 1.0;
  }
  else
  {
    /* otherwise copy SYNC-word and frame length word into output buffer*/
    ybuff[0] = xbuff[0];
    ybuff[1] = xbuff[1];

    /* ... and store data remaining data */
    for (i=2; i<lseg; i++)
      ybuff[i] = (short)(0x00FF) & (xbuff[i]);
    fer = 0.0;
  }

  /* Return the frame erasure flag */
  return(fer);
}
/* ....................... End of FER_module_stl96() ....................... */
 

/*
  ============================================================================ 

        double EID_random (unsigned long *seed);
        ~~~~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Returns a new random number, generated a linear congruential
        sequence (LCS) generator. See: Knuth, D.E. 1981: "Seminumerical
        Algorithms" vol.2 of The Art of Computer Programming; Reading,
        Mass.; Addison-Wesley.

        Parameters:    
        ~~~~~~~~~~~  
        seed: ... long seed.

        Return value:
        ~~~~~~~~~~~~~
        Returns a random number as double in the range 0..1.

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        28.Feb.92 v1.0 Release of 1st version <hf@pkinbg.uucp>
        20.Apr.92 v2.0 Change of polinomial generator to LCG
                       <hf@pkinbg.uucp>
	26.Jan.98 v3.0 Corrected bug for 64-bit operating systems (where 
                       longs have 64, not 32 bits). Implemented by
                       <simao.campos@comsat.com>, after bug reported
                       by <claude.lamblin@cnet.francetelecom.fr>

 ============================================================================
*/
double EID_random(seed)
unsigned long *seed;
{
  /* Size in bits (=size in bytes * 8) for long variables*/
   static double bits_in_long = sizeof(long) * 8;

   /* Update RNG */
   *seed = ((unsigned long)69069L * (*seed) +1L);

   /* Return random number as a double */
#ifdef WAS
   return(pow((double)2.0, (double)-32.0) * (double)(*seed));
#else
   return(pow((double)2.0, -bits_in_long) * (double)(*seed));
#endif
}
/* ....................... End of EID_random() ....................... */ 
 
 
/*
  ============================================================================ 

        long GEC_init (SCD_EID *EID, double  ber, double  gamma);
        ~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Initialise Gilbert Elliot channel model for generating bit 
        error sequences. A detailed description is given in a separate
        ITU-T document, and here only a short overview is given: 

                                       P
                  +----+   +----------->----------+   +----+
                  |    +---^-+                  +-V---+    |
              1-P |    |  G  |                  |  B  |    | 1-Q
                  |    +---^-+                  +-V---+    |
                  +--->+   +-----------<----------+   +<---+
                                       Q

        The model provides two different channel states, GOOD (G) and
        BAD (B). If the channel is in state G, the bit error rate is
        P_G, in state B the bit error rate is P_B. The channel state
        changes from G to B with probability P, and from B to G with
        probability Q. The probability, that the channel remains in the
        current state is then 1-P and 1-Q, respectively. 

        If the mean value of the bit error rate of the error pattern
        should be "ber", P and Q are computed as follows: 

                                     /      (P_B - ber) \ 
                  P   =  (1-gamma) * | 1 -  ----------- |
                                     \      (P_B - P_G) / 

                  Q   =  1  -  P  - gamma

        P_G and P_B are fixed and depend not on the selected bit error
        rate (ber). Here they are fixed as: 

                              P_G = 0.0
                              P_B = 0.5

        With the factor "gamma" the 'correlation' between adjacent bit
        errors can be defined. gamma=0 means 'no correlation' (or random
        error sequence). The more gamma runs to 1.0, the more a
        'correlation' between adjacent bit errors will occur (bursty
        error sequences). 

        Examples: 
        1) gamma = 0: (random errors) 
        
        Compute a bit error sequence of lets say 100000 bits with a bit
        error rate of 0.01. Then N auto-correlation coefficients
        (ACF(tau)) of this sequence are computed and normalized to 
        ACF(0) = 1.0. We can observe, that the remaining coefficients
        ACF(1), ACF(2), ... ACF(N) oscillate around the selected bit
        error rate of 0.01. 

        2) gamma = 0.5: (slightly bursty errors) 
        
        Do the same steps as in Example 1. Now we observe, that the
        coefficients ACF(6), ACF(7), ... ACF(N) oscillate around the
        selected bit error rate of 0.01, whereas ACF(1) ... ACF(5) 
        build a transition region. 

        The more gamma increases, the longer this transition region 
        will be! For example, with gamma=0.9 this transition ends 
        around tau=70. 

        Restrictions: 
        ------------- 
        - gamma < 0.99, otherwise no reasonable 'short term' 
                        error sequences will be produced 
        - ber  =< 0.50, otherwise erroneous sequences may be 
                        produced. 

        Parameters:    
        ~~~~~~~~~~~  
        EID: ..... pointer to EID-struct 
        ber: ..... bit error rate: 0 ... 1.0 
        gamma: ... gamma =   0,    for indepentent errors 
                            >0,<1, for blending independent burst errors 
                         =   1,    for strong burst errors 

        Return value:
        ~~~~~~~~~~~~~
        Returns a 1l on success; 0l if memory allocation failed.

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        28.Feb.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================
*/
long GEC_init(EID, ber, gamma)
SCD_EID *EID;
double  ber;
double  gamma;
{
   long i;
    double P,Q,P_G,P_B;
 

   /* simple channel model with only two states: "BAD" and "GOOD" */
    EID->nstates = 2;

  /* Allocate memory for bit error array */
    if ((EID->ber = (double *)malloc(EID->nstates*sizeof(double))) == 0L) 
      return((long)0);
 
  /* Allocate memory for the `transition matrix' */
    if ((EID->matrix = (double **)malloc(EID->nstates*(long)sizeof(double *) )) 
         == (double **)0) 
       return(0L);
 
    for (i=0; i<EID->nstates; i++)
    {
       EID->matrix[i] = (double *)malloc(EID->nstates*(long)sizeof(double));
       if (EID->matrix[i]  == (double *)0) return(0L);
    }
 
  /* Fix P_G to 0.0 and P_B to 0.5 */
    P_G = 0.0;
    P_B = 0.5;

   /* P("GOOD" --> "BAD") */
    P   = (1.0 - gamma) * 2.0 * ber;

   /* P("BAD"  --> "GOOD") */
    Q   = (1.0 - gamma) * (1.0 - 2.0 * ber);

   /* Error rate in state GOOD */
    EID->ber[0]   = P_G;

   /* Error rate in state BAD */
    EID->ber[1]   = P_B;

  /* Fill the transition matrix with normalized thresholds, which must be 
    integer values, since the random generator produces integer values from 
    0x00000000 ... 0x7FFFFFFF.  */

    EID->matrix[0][0]  = (1.0-P);               /* "GOOD" -> "GOOD" */
    EID->matrix[0][1]  =  1.0;                  /* Intervall "GOOD" -> "BAD" */
    EID->matrix[1][0]  =  Q;                    /* Intervall "BAD"  -> "GOOD" */
    EID->matrix[1][1]  =  1.0;                  /* "BAD"  -> "BAD" */
    EID->current_state = 0L;                    /* start with state "GOOD" */
 
  /* Return OK */
    return(1L);
}
/* ....................... End of GEC_init() ....................... */


/*
  ========================================================================== 

        void set_RAN_seed (SCD_EID *EID, unsigned long seed)
        ~~~~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Auxiliary function to set channel state. Modifying Channel
        Model: random generator's seed.


        Parameters:    
        ~~~~~~~~~~~  
        EID: ......... (In) Pointer to EID-struct 
        seed: ........ (In) new seed for random generator

        Return value:
        ~~~~~~~~~~~~~
        None.

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        20.Apr.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================
*/
void set_RAN_seed(EID,seed)                                                  
SCD_EID *EID;
unsigned long seed;
{
  EID->seed = seed;
}
/* ....................... End of set_RAN_seed() ....................... */


/*
  ========================================================================== 

        unsigned long get_RAN_seed (SCD_EID *EID);
        ~~~~~~~~~~~~~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Auxiliary function to support saving channel states. Get random
        number generator seed.

        Parameters:    
        ~~~~~~~~~~~  
        EID: ... pointer to EID structure with states.

        Return value:
        ~~~~~~~~~~~~~
        Returns the current seed.

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        20.Apr.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================ 
*/
unsigned long get_RAN_seed(EID)
SCD_EID *EID;
{
   return(EID->seed);
}
/* ....................... End of get_RAN_seed() ....................... */


/*
  ========================================================================== 

        void set_GEC_matrix (SCD_EID *EID, double threshold, 
        ~~~~~~~~~~~~~~~~~~~  int current_state, int next_state);        

        Description:
        ~~~~~~~~~~~~

        Auxiliary function to support setting channel states. Defines
        Transition Matrix of GEC-model.

        Parameters:    
        ~~~~~~~~~~~  
        EID: ............. (In) Pointer to EID-struct 
        threshold: ....... (In) threshold 
        current_state: ... (In) current state ('G' or 'B') 
        next_state: ...... (In) next state ('G' or 'B') 

        Return value:    none.
        ~~~~~~~~~~~~~

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        28.Feb.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================
*/
void set_GEC_matrix(EID,threshold,current_state,next_state)
SCD_EID *EID;
double threshold;
int /* char */ current_state;
int /* char */ next_state;
{
  int icurr,inext;


  if (toupper(current_state)=='G')
    icurr = 0;
  else
    icurr = 1;

  if (toupper(next_state)=='G')
    inext = 0;
  else
    inext = 1;

  EID->matrix[icurr][inext] = threshold;
}
/* ....................... End of set_GEC_matrix() ....................... */



/*
  ========================================================================== 

        double get_GEC_matrix (SCD_EID *EID, int current_state, int
        ~~~~~~~~~~~~~~~~~~~~~  next_state);
        

        Description:
        ~~~~~~~~~~~~

        Auxiliary function to support saving of channel states.
        Get the state of the Transition Matrix of GEC-model.

        Parameters:    
        ~~~~~~~~~~~  
        EID: ............. (In) Pointer to EID-struct 
        current_state: ... (In) current state ('G' or 'B') 
        next_state: ...... (In) next state ('G' or 'B') 

        Return value:
        ~~~~~~~~~~~~~
        Returns the threshold as a double.

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        20.Apr.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================
*/
double get_GEC_matrix(EID,current_state,next_state)
SCD_EID *EID;
int /* char */ current_state;
int /* char */ next_state;
{
   int icurr,inext;


   if (toupper(current_state)=='G')
     icurr = 0;
   else
     icurr = 1;

   if (toupper(next_state)=='G')
     inext = 0;
   else
     inext = 1;

   return(EID->matrix[icurr][inext]);
}
/* ....................... End of get_GEC_matrix() ....................... */


/*
  ========================================================================== 

        void set_GEC_current_state (SCD_EID *EID, int current_state);
        ~~~~~~~~~~~~~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Auxiliary function to support setting the current channel state.

        Parameters:    
        ~~~~~~~~~~~  
        EID: ............. (In) Pointer to EID-struct 
        current_state: ... (In) current state ('G' or 'B') 

        Return value:     none.
        ~~~~~~~~~~~~~

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        20.Apr.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================
*/
void set_GEC_current_state(EID,current_state)
SCD_EID *EID;
int /* char */ current_state;
{
  if (toupper(current_state)=='G')
    EID->current_state = 0L;
  else
    EID->current_state = 1L;
}
/* ................... End of set_GEC_current_state() ................... */


/*
  ========================================================================== 

        char get_GEC_current_state (SCD_EID *EID); 
        ~~~~~~~~~~~~~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Function to get the channel current state.

        Parameters:    
        ~~~~~~~~~~~  
        EID: ..... pointer to state variable;

        Return value:
        ~~~~~~~~~~~~~
        Return the current state (Good or Bad).

        Author: <hf@pkinbg.uucp>
        ~~~~~~~

        History:
        ~~~~~~~~
        20.Apr.92 v1.0 Release of 1st version <hf@pkinbg.uucp>

 ============================================================================ 
*/
char get_GEC_current_state(EID)
SCD_EID *EID; 
{
  if (EID->current_state==0L)
    return('G');
  else
    return('B');
}
/* .................... End of get_GEC_current_state() .................... */


/*
  ============================================================================ 

        BURST_EID *open_burst_eid (long index);
        ~~~~~~~~~~~~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Allocate memory for BURST_EID struct and setup the transition
        matrix according the selected bit error rate.

        Parameters:  
        ~~~~~~~~~~~
        index: ..... (In) Index for model vector parameter
			    1=0.5% 2 = 1%, 6=3%, 10=5%, 60=30% 

        Return value:
        ~~~~~~~~~~~~~
        Returns a pointer to struct BURST_EID;

        Author: 
        ~~~~~~~
	Cliff Halevi
        Bellcore, NVC 5F 319	Phone: +1-908-758-5839
        100 Schultz Drive	Fax:   +1-908-758-4502
        Red Bank, NJ 07701.	email: chh2@cc.bellcore.com

        Vijay Varma
        Bellcore NVC 3X-325	Phone: +1-908-758-2811
        331 Newman Springs Road	Fax:   +1-908-758-4371
        Red Bank, NJ 07701.	email: varma@nyquist.bellcore.com

	Gerhard Schroeder .....

        History:
        ~~~~~~~~
        04/Aug/93 v.1.0 Created

 ============================================================================
*/


BURST_EID   *open_burst_eid (index)
long index;
{
  BURST_EID *burst_eid;
  time_t t1;
  long i;

  /* Allocate EID structure */
   if ((burst_eid=(BURST_EID *)malloc(sizeof(BURST_EID)))==0L) 
	    return((BURST_EID *)0);

  /* Preset of the random generator seed with current system time */
#ifdef PORT_TEST
    t1 = 314159265;
#else
    time(&t1);
#endif

    burst_eid->seedptr = (unsigned long) t1;
    for (i=0;  i<MODEL_SIZE;  i++)
    {
      burst_eid->internal[i] = 0;
    }
    burst_eid->s_new = 0;
    
    burst_eid->index = index;
	if((burst_eid->index>60) || (burst_eid->index < 1)){
	    fprintf(stderr, " Using default error rate of 1 %%\n");
	    burst_eid->index = 2;
    }

    prob[0]=bfer_comp(burst_eid->index);

#ifdef DEBUG 
    /*	 
    **  Print out the selected coefficient vector
    */        
    for (i=0;  i< MODEL_SIZE;  i++){
      printf ("p[%2d]=%f\n", i,prob[i]);
    }
#endif

   /* Return BURST_EID structure */
    return (burst_eid);
}
/* .................... End of open_burst_eid() .................... */


/*
  ============================================================================ 

        double FER_generator_burst (BURST_EID *state);
        ~~~~~~~~~~~~~~~~~~~~~~~~~~

        Description:
        ~~~~~~~~~~~~

        Generate a double indicating whether the current frame should be
        considered good or bad..

        Parameters:  
        ~~~~~~~~~~~
        state: ..... (In/Out) State variable describing the channel.

        Return value:
        ~~~~~~~~~~~~~
        Returns 0 on a good-frame indication and 1 if the current frame
        is to be considered a bad frame.

        Author: 
        ~~~~~~~
	Gerhard Schroeder
	FI-DBP Telekom            Tel: 	+49-6151-83-3973
	Am Kavalleriesand 3       Fax: 	+49-6151-89-5234
	D-64295 Darmstadt         E-mail: gerhard.schroeder@fz13.fz.dbp.de
	GERMANY

        History:
        ~~~~~~~~
        04/Aug/93 v1.0 Created
	26.Jan.98 v3.0 Corrected bug for 64-bit operating systems (where 
                       longs have 64, not 32 bits). Implemented by
                       <simao.campos@comsat.com>, after bug reported
                       by <claude.lamblin@cnet.francetelecom.fr>

 ============================================================================
*/
double FER_generator_burst (state)
BURST_EID *state;
{
  long    in;
  double  ran, aux; 
  /* Size in bits (=size in bytes * 8) for long variables*/
  static double bits_in_long = sizeof(long) * 8;

  /* UGST Random number generator */
  state->seedptr = ((unsigned long)69069L * (state->seedptr) +1L);

  /* Return random number as a double in range [0...1] */
#ifdef WAS
  ran = (pow((double)2.0, (double)-32.0) * (double)(state->seedptr));
#else
  ran = (pow((double)2.0, -bits_in_long) * (double)(state->seedptr));
#endif
   aux = prob[state->s_new];
  in = floor(ran + aux);	/* in = 0 indicates good frame */
                                /* in = 1 indicates bad frame */
  if (in == 0)                    /* If frame is good, */
  {
    /* increment frame counter for state */
    state->internal[state->s_new]++;
    if (state->s_new != 0) state->internal[0]++;
    state->s_new = 0;
  }
  else
    state->s_new++;      /* Else, if frame is bad, increment state */

  return((double) in);
}
/* .................... End of FER_generator_burst() .................... */


/* ********************************************************************* */
/* ********************************************************************* */
/* ************************* ADDED AFTER STL96 ************************* */
/* ********************************************************************* */
/* ********************************************************************* */

/* 
   -------------------------------------------------------------------------
   void reset_burst_eid (BURST_EID *burst_eid);
   ~~~~~~~~~~~~~~~~~~~~

   Reset a BURST_EID structure. By default, only counters are reset;
   if the symbol RESET_SEED_AS_WELL is defined at compilation time,
   the seed is also reset. However, this is not recommended.

   Parameter:
   ~~~~~~~~~~
   burst_eid: .... BURST_EID structure (previously initialized by 
                   open_burst_eid) to be reset.

   Return value: None.
   ~~~~~~~~~~~~~

   Original author: <simao.campos@comsat.com>
   ~~~~~~~~~~~~~~~~

   History:
   ~~~~~~~~
   15.Aug.97  v.1.0  Created.
   -------------------------------------------------------------------------
 */
BURST_EID *reset_burst_eid(burst_eid)
BURST_EID *burst_eid;
{
#ifdef RESET_SEED_AS_WELL
  time_t t1;
#endif
  long i;

  /* If state variable was not initialized, return a null pointer */
  if (burst_eid == (BURST_EID *)0)
    return((BURST_EID *)0);

#ifdef RESET_SEED_AS_WELL
  /* Reset random generator seed */
# ifdef PORT_TEST
  t1 = 314159265;
# else
  time(&t1);
# endif
  burst_eid->seedptr = (unsigned long) t1;
#endif  /* RESET_SEED_AS_WELL */

  /* Reset internal counters */
  for (i=0;  i<MODEL_SIZE;  i++)
  {
    burst_eid->internal[i] = 0;
  }
  burst_eid->s_new = 0;

  /* Return pointer to BURST_EID state variable */
  return (burst_eid);
}
/* ...................... End of reset_burst_eid() ...................... */

/* bfer_comp for Bellcore model
   compute steady state value of P0 in state transition prob matrix,
   given that P1..P9 are fixed 
   input index, gives desired total output FER(Pe) */ 

double bfer_comp(long index)
{
	double Pe;
	double num,den,prod;
	int i,k;

	Pe = (double)index/200.0; /* Pe quantized in 0.5% steps due to legacy state saving method */

	num=(1/(1-Pe))-1.0;
    den=0.0;
	for (i=1;i<=(MODEL_SIZE-2);i++){
	   prod=1.0;
	   for(k=1;k<=i;k++){
         prod *= prob[k];
	   }
	   den+=prod;
	 }
	 den=den+1.0;

 return (num/den); 
}
